module IrrigationMod

  !-----------------------------------------------------------------------
  ! !DESCRIPTION:
  ! Calculates irrigation flux.
  !
  ! Usage:
  !
  !   - Call CalcIrrigationNeeded in order to compute whether and how much irrigation is
  !     needed for the next call to CalcIrrigationFluxes. This should be called once per
  !     timestep.
  ! 
  !   - Call CalcIrrigationFluxes in order to calculate irrigation withdrawal and
  !     application fluxes. It is acceptable for this to be called earlier in the timestep
  !     than CalcIrrigationNeeded.
  !
  ! General design notes:
  !
  !   In principle, CalcIrrigationFluxes and CalcIrrigationNeeded could be combined into a
  !   single routine. Their separation is largely for historical reasons: In the past,
  !   irrigation depended on btran, and qflx_irrig is needed earlier in the driver loop
  !   than when btran becomes available. (And qflx_irrig is also used late in the driver
  !   loop - so it wouldn't work, for example, to calculate qflx_irrig after btran is
  !   computed, and then save it on the restart file for the next iteration of the driver
  !   loop: then the uses of qflx_irrig early and late in the driver loop would be
  !   inconsistent.)
  !
  !   Now that we no longer have a dependency on btran, we could call CalcIrrigationNeeded
  !   before the first time qflx_irrig is needed. Thus, there might be some advantage to
  !   combining CalcIrrigationFluxes and CalcIrrigationNeeded - or at least calling these
  !   two routines from the same place.  In particular: this separation of the irrigation
  !   calculation into two routines that are done at different times in the driver loop
  !   makes it harder and less desirable to nest the irrigation object within some other
  !   object: Doing so might make it harder to do the two separate steps at the right
  !   time, and would lead to less clarity about how these two steps are ordered with
  !   respect to the rest of the driver loop. So if we start trying to create a hierarchy
  !   of objects in CLM, we may want to rework this design. And we may want to rework it
  !   to keep things simpler anyway, even if we aren't nesting objects. Note that this
  !   rework would change behavior slightly, because irrigation would be applied in the
  !   same time step that CalcIrrigationNeeded first determines it's needed, rather than
  !   waiting until the following time step.
  !
  ! Design notes regarding where we calculate supply-based irrigation limitations:
  !
  !   We compute supply-based irrigation limitations differently for rivers (volr) vs.
  !   groundwater.
  !
  !   For volr-based limitation, we apply the limit in CalcIrrigationNeeded rather than in
  !   CalcIrrigationFluxes. This is to avoid problems that arise due to evolving volr -
  !   e.g., if we had calculated irrigation demand such that volr was just barely big
  !   enough, then we apply irrigation for some time steps, then we get an updated volr
  !   which is lower due to these irrigation withdrawals - if we applied the volr
  !   limitation in CalcIrrigationFluxes, this would impose an unnecessary limit on the
  !   irrigation to be applied for the rest of this day. It's hard to see a clean way
  !   around this problem, given that volr is updated in some but not all CLM time
  !   steps. So we impose the limit on irrigation deficit, rather than on the
  !   time step-by-time-step flux. The downside here is that it's somewhat more likely
  !   that irrigation would try to draw volr negative, if volr has been evolving for other
  !   reasons.
  !
  !   For groundwater-based limitation, in contrast, we can rely on having an up-to-date
  !   view of available groundwater, so the argument against having the volr limitation
  !   in CalcIrrigationFluxes doesn't apply, and that seems like a more robust place to
  !   do this limitation.
  !
  !   The key difference is that, for volr, we don't have an easy way to determine if this
  !   is a time step when we have an updated volr value, or if we still have volr from a
  !   few time steps ago (since ROF isn't coupled in every time step). So in that case,
  !   having the possibility that we'd exceed available volr seems like the lesser of
  !   evils.
  !
  ! !USES:
#include "shr_assert.h"
  use shr_kind_mod     , only : r8 => shr_kind_r8
  use decompMod        , only : bounds_type, get_proc_global
  use decompMod        , only : subgrid_level_gridcell, subgrid_level_column, subgrid_level_patch
  use shr_log_mod      , only : errMsg => shr_log_errMsg
  use abortutils       , only : endrun
  use clm_instur       , only : irrig_method
  use pftconMod        , only : pftcon
  use clm_varctl       , only : iulog
  use clm_varcon       , only : isecspday, denh2o, spval, ispval
  use clm_varpar       , only : nlevsoi, nlevgrnd
  use clm_time_manager , only : get_step_size
  use SoilHydrologyMod , only : CalcIrrigWithdrawals
  use SoilHydrologyType, only : soilhydrology_type
  use SoilStateType    , only : soilstate_type
  use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
  use WaterType        , only : water_type
  use WaterFluxBulkType, only : waterfluxbulk_type
  use WaterFluxType    , only : waterflux_type
  use WaterStateBulkType, only : waterstatebulk_type
  use WaterStateType   , only : waterstate_type
  use WaterTracerUtils , only : CalcTracerFromBulk, CalcTracerFromBulkFixedRatio
  use GridcellType     , only : grc
  use ColumnType       , only : col                
  use PatchType        , only : patch                
  use subgridAveMod    , only : p2c, c2g
  use filterColMod     , only : filter_col_type, col_filter_from_logical_array
  !
  implicit none
  private

  ! !PUBLIC TYPES:
  
  ! This type is public (and its components are public, too) to aid unit testing
  type, public :: irrigation_params_type
     ! Minimum LAI for irrigation
     real(r8) :: irrig_min_lai

     ! Time of day to check whether we need irrigation, seconds (0 = midnight). 
     ! We start applying the irrigation in the time step FOLLOWING this time, 
     ! since we won't begin irrigating until the next call to CalcIrrigationFluxes
     integer  :: irrig_start_time

     ! Desired amount of time to irrigate per day (sec). Actual time may 
     ! differ if this is not a multiple of dtime. Irrigation won't work properly 
     ! if dtime > secsperday
     integer  :: irrig_length

     ! Target soil matric potential for irrigation (mm)
     !
     ! When we irrigate, we aim to bring the total soil moisture in the top (irrig_depth)
     ! m of soil up to this level.
     !
     ! (Note: When we convert this to a relative saturation, we ensure that the relative
     ! saturation target is bounded by [0,1].)
     real(r8) :: irrig_target_smp

     ! Soil depth to which we measure for irrigation (m)
     real(r8) :: irrig_depth

     ! Determines soil moisture threshold at which we irrigate. If
     ! h2osoi_liq_wilting_point is the soil moisture level at wilting point and
     ! h2osoi_liq_target is the soil moisture level at the target irrigation level (given
     ! by irrig_target_smp), then the threshold at which we irrigate is
     !     h2osoi_liq_wilting_point +
     !          irrig_threshold_fraction*(h2osoi_liq_target - h2osoi_liq_wilting_point)
     ! A value of 1 means that we irrigate whenever soil moisture falls below the target
     ! A value of 0 means that we only irrigate when soil moisture falls below the
     ! wilting point
     real(r8) :: irrig_threshold_fraction

     ! Threshold for river water volume below which irrigation is shut off, if
     ! limit_irrigation is .true. (fraction of available river water). A threshold of 0
     ! means allow all river water to be used; a threshold of 0.1 means allow 90% of the
     ! river volume to be used; etc.
     real(r8) :: irrig_river_volume_threshold

     ! Whether irrigation is limited based on river storage. This only applies if ROF is
     ! enabled (i.e., rof_prognostic is .true.) - otherwise we don't limit irrigation,
     ! regardless of the value of this flag.
     logical :: limit_irrigation_if_rof_enabled

     ! use groundwater supply for irrigation (in addition to surface water)
     logical :: use_groundwater_irrigation

     ! Irrigation method to use if not set on the surface dataset
     ! (This won't be needed once we can rely on irrig_method always being on the surface
     ! dataset, but it is useful until then - particularly for testing. See also
     ! https://github.com/ESCOMP/ctsm/issues/565.)
     integer :: irrig_method_default

  end type irrigation_params_type

  type, public :: irrigation_type
     private

     ! Private data members; set in initialization:
     type(irrigation_params_type) :: params
     integer :: dtime                ! land model time step (sec)
     integer :: irrig_nsteps_per_day ! number of time steps per day in which we irrigate
     real(r8), pointer :: relsat_wilting_point_col(:,:) ! relative saturation at which smp = wilting point [col, nlevsoi]
     real(r8), pointer :: relsat_target_col(:,:)        ! relative saturation at which smp is at the irrigation target [col, nlevsoi]
     integer , pointer :: irrig_method_patch          (:) ! patch irrigation application method

     ! Private data members; time-varying:
     real(r8), pointer :: sfc_irrig_rate_patch        (:) ! current irrigation rate from surface water [mm/s]
     real(r8), pointer :: irrig_rate_demand_patch     (:) ! current irrigation rate, neglecting surface water source limitation [mm/s]
     integer , pointer :: n_irrig_steps_left_patch    (:) ! number of time steps for which we still need to irrigate today (if 0, ignore)
     real(r8), pointer :: qflx_irrig_demand_patch     (:) ! irrigation flux neglecting surface water source limitation [mm/s]

   contains
     ! Public routines
     ! COMPILER_BUG(wjs, 2014-10-15, pgi 14.7) Add an "Irrigation" prefix to some  generic routines like "Init"
     ! (without this workaround, pgi compilation fails in restFileMod)
     procedure, public :: Init => IrrigationInit
     procedure, public :: Restart
     procedure, public :: CalcIrrigationFluxes
     procedure, public :: CalcIrrigationNeeded
     procedure, public :: UseGroundwaterIrrigation ! Returns true if groundwater irrigation enabled
     procedure, public :: Clean => IrrigationClean ! deallocate memory

     ! Public simply to support unit testing; should not be used from CLM code
     procedure, public :: InitForTesting ! version of Init meant for unit testing
     procedure, public :: RelsatToH2osoi ! convert from relative saturation to kg/m2 water for a single column and layer
     procedure, public :: WrapCalcIrrigWithdrawals ! Wrap the call to CalcIrrigWithdrawals

     ! Private routines
     procedure, private :: ReadNamelist
     procedure, private :: CheckNamelistValidity   ! Check for validity of input parameters
     procedure, private :: InitAllocate => IrrigationInitAllocate
     procedure, private :: InitHistory => IrrigationInitHistory
     procedure, private :: InitCold => IrrigationInitCold
     procedure, private :: CalcBulkWithdrawals      ! calculate irrigation withdrawals for bulk water
     procedure, private :: CalcOneTracerWithdrawals ! calculate irrigation withdrawals for one water tracer
     procedure, private :: CalcTotalGWUnconIrrig    ! calculate total irrigation withdrawal flux from the unconfined aquifer, for either bulk or one water tracer
     procedure, private :: CalcApplicationFluxes    ! calculate irrigation application fluxes for either bulk water or a single tracer
     procedure, private :: CalcIrrigNstepsPerDay    ! given dtime, calculate irrig_nsteps_per_day
     procedure, private :: SetIrrigMethod           ! set irrig_method_patch based on surface dataset
     procedure, private :: PointNeedsCheckForIrrig  ! whether a given point needs to be checked for irrigation now
     procedure, private :: CalcDeficitVolrLimited   ! calculate deficit limited by river volume for each patch
  end type irrigation_type

  interface irrigation_params_type
     module procedure irrigation_params_constructor
  end interface irrigation_params_type

  ! Soil matric potential at wilting point (mm)
  !
  ! There is no reason to make this a tunable parameter, because the behavior it governs
  ! (the trigger for irrigation) can be tuned via other parameters.
  !
  ! TODO(wjs, 2016-09-08) It looks like there is other code in CLM that also uses an
  ! assumed wilting point (CNRootDynMod, maybe others). We should probably make this a
  ! shared parameter, e.g., in clm_varcon.
  real(r8), parameter, private :: wilting_point_smp = -150000._r8

  ! Conversion factors
  real(r8), parameter :: m3_over_km2_to_mm = 1.e-3_r8
     
  ! Irrigation methods
  integer, parameter, public  :: irrig_method_unset = 0
  ! Drip is defined here as irrigation applied directly to soil surface
  integer, parameter, public :: irrig_method_drip = 1
  ! Sprinkler is applied directly to canopy
  integer, parameter, public :: irrig_method_sprinkler = 2

  character(len=*), parameter, private :: sourcefile = &
       __FILE__
  
contains

  ! ========================================================================
  ! Constructors
  ! ========================================================================

  !-----------------------------------------------------------------------
  function irrigation_params_constructor(irrig_min_lai, &
       irrig_start_time, irrig_length, &
       irrig_target_smp, &
       irrig_depth, irrig_threshold_fraction, irrig_river_volume_threshold, &
       limit_irrigation_if_rof_enabled, use_groundwater_irrigation, &
       irrig_method_default) &
       result(this)
    !
    ! !DESCRIPTION:
    ! Create an irrigation_params instance
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    type(irrigation_params_type) :: this  ! function result
    real(r8), intent(in) :: irrig_min_lai
    integer , intent(in) :: irrig_start_time
    integer , intent(in) :: irrig_length
    real(r8), intent(in) :: irrig_target_smp
    real(r8), intent(in) :: irrig_depth
    real(r8), intent(in) :: irrig_threshold_fraction
    real(r8), intent(in) :: irrig_river_volume_threshold
    logical , intent(in) :: limit_irrigation_if_rof_enabled
    logical , intent(in) :: use_groundwater_irrigation
    integer , intent(in) :: irrig_method_default
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'irrigation_params_constructor'
    !-----------------------------------------------------------------------
    
    this%irrig_min_lai = irrig_min_lai
    this%irrig_start_time = irrig_start_time
    this%irrig_length = irrig_length
    this%irrig_depth = irrig_depth
    this%irrig_target_smp = irrig_target_smp
    this%irrig_threshold_fraction = irrig_threshold_fraction
    this%irrig_river_volume_threshold = irrig_river_volume_threshold
    this%limit_irrigation_if_rof_enabled = limit_irrigation_if_rof_enabled
    this%use_groundwater_irrigation = use_groundwater_irrigation
    this%irrig_method_default = irrig_method_default

  end function irrigation_params_constructor


  ! ========================================================================
  ! Infrastructure routines (initialization, restart, etc.)
  ! ========================================================================
  
  !------------------------------------------------------------------------
  subroutine IrrigationInit(this, bounds, NLFilename, &
       soilstate_inst, soil_water_retention_curve, &
       use_aquifer_layer)
    use SoilStateType , only : soilstate_type

    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    character(len=*)       , intent(in)    :: NLFilename ! Namelist filename
    type(soilstate_type)   , intent(in)    :: soilstate_inst
    class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve
    logical                , intent(in)    :: use_aquifer_layer ! whether an aquifer layer is used in this run

    call this%ReadNamelist(NLFilename, use_aquifer_layer)
    call this%InitAllocate(bounds)
    call this%InitHistory(bounds)
    call this%InitCold(bounds, soilstate_inst, soil_water_retention_curve)
  end subroutine IrrigationInit

  !-----------------------------------------------------------------------
  subroutine InitForTesting(this, bounds, params, dtime, &
       relsat_wilting_point, relsat_target)
    !
    ! !DESCRIPTION:
    ! Does initialization needed for unit testing. Allows caller to prescribe values of
    ! some internal variables.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(irrigation_type)       , intent(inout) :: this
    type(bounds_type)            , intent(in)    :: bounds
    type(irrigation_params_type) , intent(in)    :: params
    integer                      , intent(in)    :: dtime ! model time step (sec)
    real(r8)                     , intent(in)    :: relsat_wilting_point( bounds%begc: , 1: ) ! relative saturation at which smp = irrig_wilting_point_smp [col, nlevsoi]
    real(r8)                     , intent(in)    :: relsat_target( bounds%begc: , 1: ) ! relative saturation at which smp = irrig_target_smp [col, nlevsoi]
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'InitForTesting'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(relsat_wilting_point) == (/bounds%endc, nlevsoi/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(relsat_target) == (/bounds%endc, nlevsoi/)), sourcefile, __LINE__)

    call this%InitAllocate(bounds)
    this%params = params
    this%dtime = dtime
    call this%SetIrrigMethod(bounds)
    this%irrig_nsteps_per_day = this%CalcIrrigNstepsPerDay(dtime)
    this%relsat_wilting_point_col(:,:) = relsat_wilting_point(:,:)
    this%relsat_target_col(:,:) = relsat_target(:,:)

  end subroutine InitForTesting

  !-----------------------------------------------------------------------
  subroutine ReadNamelist(this, NLFilename, use_aquifer_layer)
    !
    ! !DESCRIPTION:
    ! Read the irrigation namelist
    !
    ! !USES:
    use fileutils      , only : getavu, relavu, opnfil
    use shr_nl_mod     , only : shr_nl_find_group_name
    use spmdMod        , only : masterproc, mpicom
    use shr_mpi_mod    , only : shr_mpi_bcast
    use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    character(len=*), intent(in) :: NLFilename ! Namelist filename
    logical, intent(in) :: use_aquifer_layer    ! whether an aquifer layer is used in this run
    !
    ! !LOCAL VARIABLES:

    ! temporary variables corresponding to the components of irrigation_params_type
    real(r8) :: irrig_min_lai
    integer  :: irrig_start_time
    integer  :: irrig_length
    real(r8) :: irrig_target_smp
    real(r8) :: irrig_depth
    real(r8) :: irrig_threshold_fraction
    real(r8) :: irrig_river_volume_threshold
    logical  :: limit_irrigation_if_rof_enabled
    logical  :: use_groundwater_irrigation
    character(len=64) :: irrig_method_default
    integer  :: irrig_method_default_int

    integer :: ierr                 ! error code
    integer :: unitn                ! unit for namelist file
    character(len=*), parameter :: nmlname = 'irrigation_inparm'

    character(len=*), parameter :: subname = 'ReadNamelist'
    !-----------------------------------------------------------------------

    namelist /irrigation_inparm/ irrig_min_lai, irrig_start_time, irrig_length, &
         irrig_target_smp, irrig_depth, irrig_threshold_fraction, &
         irrig_river_volume_threshold, limit_irrigation_if_rof_enabled, &
         use_groundwater_irrigation, irrig_method_default
    
    ! Initialize options to garbage defaults, forcing all to be specified explicitly in
    ! order to get reasonable results
    irrig_min_lai = nan
    irrig_start_time = 0
    irrig_length = 0
    irrig_target_smp = nan
    irrig_depth = nan
    irrig_threshold_fraction = nan
    irrig_river_volume_threshold = nan
    limit_irrigation_if_rof_enabled = .false.
    use_groundwater_irrigation = .false.
    irrig_method_default = ' '

    if (masterproc) then
       unitn = getavu()
       write(iulog,*) 'Read in '//nmlname//'  namelist'
       call opnfil (NLFilename, unitn, 'F')
       call shr_nl_find_group_name(unitn, nmlname, status=ierr)
       if (ierr == 0) then
          read(unitn, nml=irrigation_inparm, iostat=ierr)
          if (ierr /= 0) then
             call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
          end if
       else
          call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__))
       end if
       call relavu( unitn )
    end if

    call shr_mpi_bcast(irrig_min_lai, mpicom)
    call shr_mpi_bcast(irrig_start_time, mpicom)
    call shr_mpi_bcast(irrig_length, mpicom)
    call shr_mpi_bcast(irrig_target_smp, mpicom)
    call shr_mpi_bcast(irrig_depth, mpicom)
    call shr_mpi_bcast(irrig_threshold_fraction, mpicom)
    call shr_mpi_bcast(irrig_river_volume_threshold, mpicom)
    call shr_mpi_bcast(limit_irrigation_if_rof_enabled, mpicom)
    call shr_mpi_bcast(use_groundwater_irrigation, mpicom)
    call shr_mpi_bcast(irrig_method_default, mpicom)

    call translate_irrig_method_default

    this%params = irrigation_params_type( &
         irrig_min_lai = irrig_min_lai, &
         irrig_start_time = irrig_start_time, &
         irrig_length = irrig_length, &
         irrig_target_smp = irrig_target_smp, &
         irrig_depth = irrig_depth, &
         irrig_threshold_fraction = irrig_threshold_fraction, &
         irrig_river_volume_threshold = irrig_river_volume_threshold, &
         limit_irrigation_if_rof_enabled = limit_irrigation_if_rof_enabled, &
         use_groundwater_irrigation = use_groundwater_irrigation, &
         irrig_method_default = irrig_method_default_int)

    if (masterproc) then
       write(iulog,*) ' '
       write(iulog,*) nmlname//' settings:'
       ! Write settings one-by-one rather than with a nml write because
       ! irrig_river_volume_threshold may be NaN
       write(iulog,*) 'irrig_min_lai = ', irrig_min_lai
       write(iulog,*) 'irrig_start_time = ', irrig_start_time
       write(iulog,*) 'irrig_length = ', irrig_length
       write(iulog,*) 'irrig_target_smp = ', irrig_target_smp
       write(iulog,*) 'irrig_depth = ', irrig_depth
       write(iulog,*) 'irrig_threshold_fraction = ', irrig_threshold_fraction
       write(iulog,*) 'limit_irrigation_if_rof_enabled = ', limit_irrigation_if_rof_enabled
       if (limit_irrigation_if_rof_enabled) then
          write(iulog,*) 'irrig_river_volume_threshold = ', irrig_river_volume_threshold
       end if
       write(iulog,*) 'use_groundwater_irrigation = ', use_groundwater_irrigation
       write(iulog,*) 'irrig_method_default = ', irrig_method_default
       write(iulog,*) ' '

       call this%CheckNamelistValidity(use_aquifer_layer)
    end if

  contains
    subroutine translate_irrig_method_default
      select case (irrig_method_default)
      case ('drip')
         irrig_method_default_int = irrig_method_drip
      case ('sprinkler')
         irrig_method_default_int = irrig_method_sprinkler
      case default
         write(iulog,*) 'ERROR: unknown irrig_method_default: ', trim(irrig_method_default)
         call endrun('Unknown irrig_method_default')
      end select
    end subroutine translate_irrig_method_default

  end subroutine ReadNamelist

  !-----------------------------------------------------------------------
  subroutine CheckNamelistValidity(this, use_aquifer_layer)
    !
    ! !DESCRIPTION:
    ! Check for validity of input parameters.
    !
    ! Assumes that the inputs have already been packed into 'this%params'.
    !
    ! Only needs to be called by the master task, since parameters are the same for all
    ! tasks.
    !
    ! !ARGUMENTS:
    class(irrigation_type), intent(in) :: this
    logical, intent(in) :: use_aquifer_layer    ! whether an aquifer layer is used in this run
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'CheckNamelistValidity'
    !-----------------------------------------------------------------------

    associate( &
         irrig_min_lai => this%params%irrig_min_lai, &
         irrig_start_time => this%params%irrig_start_time, &
         irrig_length => this%params%irrig_length, &
         irrig_target_smp => this%params%irrig_target_smp, &
         irrig_depth => this%params%irrig_depth, &
         irrig_threshold_fraction => this%params%irrig_threshold_fraction, &
         irrig_river_volume_threshold => this%params%irrig_river_volume_threshold, &
         use_groundwater_irrigation => this%params%use_groundwater_irrigation, &
         limit_irrigation_if_rof_enabled => this%params%limit_irrigation_if_rof_enabled)

    if (irrig_min_lai < 0._r8) then
       write(iulog,*) ' ERROR: irrig_min_lai must be >= 0'
       write(iulog,*) ' irrig_min_lai = ', irrig_min_lai
       call endrun(msg=' ERROR: irrig_min_lai must be >= 0 ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_start_time < 0 .or. irrig_start_time >= isecspday) then
       write(iulog,*) ' ERROR: irrig_start_time must be >= 0 and < ', isecspday
       write(iulog,*) ' irrig_start_time = ', irrig_start_time
       call endrun(msg=' ERROR: irrig_start_time out of bounds ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_length <= 0 .or. irrig_length > isecspday) then
       write(iulog,*) ' ERROR: irrig_length must be > 0 and <= ', isecspday
       write(iulog,*) ' irrig_length = ', irrig_length
       call endrun(msg=' ERROR: irrig_length out of bounds ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_target_smp >= 0._r8) then
       write(iulog,*) ' ERROR: irrig_target_smp must be negative'
       write(iulog,*) ' irrig_target_smp = ', irrig_target_smp
       call endrun(msg=' ERROR: irrig_target_smp must be negative ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_target_smp < wilting_point_smp) then
       write(iulog,*) ' ERROR: irrig_target_smp must be >= wilting_point_smp'
       write(iulog,*) ' irrig_target_smp (from namelist) = ', irrig_target_smp
       write(iulog,*) ' wilting_point_smp (hard-coded) = ', wilting_point_smp
       call endrun(msg=' ERROR: irrig_target_smp must be >= wilting_point_smp ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_depth < 0._r8) then
       write(iulog,*) ' ERROR: irrig_depth must be > 0'
       write(iulog,*) ' irrig_depth = ', irrig_depth
       call endrun(msg=' ERROR: irrig_depth must be > 0 ' // errMsg(sourcefile, __LINE__))
    end if

    if (irrig_threshold_fraction < 0._r8 .or. irrig_threshold_fraction > 1._r8) then
       write(iulog,*) ' ERROR: irrig_threshold_fraction must be between 0 and 1'
       write(iulog,*) ' irrig_threshold_fraction = ', irrig_threshold_fraction
       call endrun(msg=' ERROR: irrig_threshold_fraction must be between 0 and 1 ' // &
            errMsg(sourcefile, __LINE__))
    end if

    if (limit_irrigation_if_rof_enabled) then
       if (irrig_river_volume_threshold < 0._r8 .or. irrig_river_volume_threshold > 1._r8) then
          write(iulog,*) ' ERROR: irrig_river_volume_threshold must be between 0 and 1'
          write(iulog,*) ' irrig_river_volume_threshold = ', irrig_river_volume_threshold
          call endrun(msg=' ERROR: irrig_river_volume_threshold must be between 0 and 1 ' // &
               errMsg(sourcefile, __LINE__))
       end if
    end if

    if (use_groundwater_irrigation .and. .not. limit_irrigation_if_rof_enabled) then
       write(iulog,*) ' ERROR: use_groundwater_irrigation only makes sense if limit_irrigation_if_rof_enabled is set.'
       write(iulog,*) '(If limit_irrigation_if_rof_enabled is .false., then groundwater extraction will never be invoked.)'
       call endrun(msg=' ERROR: use_groundwater_irrigation only makes sense if limit_irrigation_if_rof_enabled is set' // &
            errMsg(sourcefile, __LINE__))
    end if

    if (use_aquifer_layer .and. use_groundwater_irrigation) then
          write(iulog,*) ' ERROR: use_groundwater_irrigation and use_aquifer_layer may not be used simultaneously'
          call endrun(msg=' ERROR: use_groundwater_irrigation and use_aquifer_layer cannot both be set to true' // &
               errMsg(sourcefile, __LINE__))
    end if

    end associate

  end subroutine CheckNamelistValidity


  !-----------------------------------------------------------------------
  subroutine IrrigationInitAllocate(this, bounds)
    !
    ! !DESCRIPTION:
    ! Initialize irrigation data structure
    !
    ! !USES:
    use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    !
    ! !LOCAL VARIABLES:
    integer :: begp, endp
    integer :: begc, endc

    character(len=*), parameter :: subname = 'InitAllocate'
    !-----------------------------------------------------------------------

    begp = bounds%begp; endp= bounds%endp
    begc = bounds%begc; endc= bounds%endc

    allocate(this%qflx_irrig_demand_patch (begp:endp))              ; this%qflx_irrig_demand_patch  (:)   = nan
    allocate(this%relsat_wilting_point_col    (begc:endc,nlevsoi))  ; this%relsat_wilting_point_col     (:,:) = nan
    allocate(this%relsat_target_col           (begc:endc,nlevsoi))  ; this%relsat_target_col            (:,:) = nan
    allocate(this%sfc_irrig_rate_patch        (begp:endp))          ; this%sfc_irrig_rate_patch         (:)   = nan
    allocate(this%irrig_rate_demand_patch     (begp:endp))          ; this%irrig_rate_demand_patch      (:)   = nan
    allocate(this%irrig_method_patch          (begp:endp))          ; this%irrig_method_patch           (:)   = ispval
    allocate(this%n_irrig_steps_left_patch    (begp:endp))          ; this%n_irrig_steps_left_patch     (:)   = 0

  end subroutine IrrigationInitAllocate

  !-----------------------------------------------------------------------
  subroutine IrrigationInitHistory(this, bounds)
    !
    ! !DESCRIPTION:
    ! Initialize irrigation history fields
    !
    ! !USES:
    use histFileMod  , only : hist_addfld1d
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    !
    ! !LOCAL VARIABLES:
    integer :: begp, endp
    
    character(len=*), parameter :: subname = 'InitHistory'
    !-----------------------------------------------------------------------

    begp = bounds%begp; endp= bounds%endp

    this%qflx_irrig_demand_patch(begp:endp) = spval
    call hist_addfld1d (fname='QIRRIG_DEMAND', units='mm/s', &
         avgflag='A', long_name='irrigation demand', &
         ptr_patch=this%qflx_irrig_demand_patch, default='inactive')

  end subroutine IrrigationInitHistory

  !-----------------------------------------------------------------------
  subroutine IrrigationInitCold(this, bounds, soilstate_inst, soil_water_retention_curve)
    !
    ! !DESCRIPTION:
    ! Do cold-start initialization for irrigation data structure
    !
    ! !USES:
    use SoilStateType , only : soilstate_type
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    type(soilstate_type)   , intent(in)    :: soilstate_inst
    class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve
    !
    ! !LOCAL VARIABLES:
    integer :: c ! col index
    integer :: j ! level index

    character(len=*), parameter :: subname = 'InitCold'
    !-----------------------------------------------------------------------

    do j = 1, nlevsoi
       do c = bounds%begc, bounds%endc
          call soil_water_retention_curve%soil_suction_inverse( &
               c = c, &
               j = j, &
               smp_target = wilting_point_smp, &
               soilstate_inst = soilstate_inst, &
               s_target = this%relsat_wilting_point_col(c,j))

          call soil_water_retention_curve%soil_suction_inverse( &
               c = c, &
               j = j, &
               smp_target = this%params%irrig_target_smp, &
               soilstate_inst = soilstate_inst, &
               s_target = this%relsat_target_col(c,j))

          ! Make sure relative saturation targets are bounded by [0,1]
          !
          ! NOTE(wjs, 2016-11-17) These targets can be > 1 if smp_target is too small of
          ! a negative value; we want to force the target to a relative saturation of 1
          ! in that case. I don't see how these targets could end up negative, though I
          ! had a note that I saw negative values at one point. In practice, for
          ! reasonable parameter values, it seems that these min and max functions have
          ! no effect.
          this%relsat_wilting_point_col(c,j) = min(this%relsat_wilting_point_col(c,j), 1._r8)
          this%relsat_wilting_point_col(c,j) = max(this%relsat_wilting_point_col(c,j), 0._r8)
          this%relsat_target_col(c,j) = min(this%relsat_target_col(c,j), 1._r8)
          this%relsat_target_col(c,j) = max(this%relsat_target_col(c,j), 0._r8)
       end do
    end do

    call this%SetIrrigMethod(bounds)
       
    this%dtime = get_step_size()
    this%irrig_nsteps_per_day = this%CalcIrrigNstepsPerDay(this%dtime)

    this%qflx_irrig_demand_patch(bounds%begp:bounds%endp) = 0._r8
    
  end subroutine IrrigationInitCold

  !-----------------------------------------------------------------------
  pure function CalcIrrigNstepsPerDay(this, dtime) result(irrig_nsteps_per_day)
    !
    ! !DESCRIPTION:
    ! Given dtime (sec), determine number of irrigation steps per day
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    integer :: irrig_nsteps_per_day  ! function result
    class(irrigation_type) , intent(in) :: this
    integer                , intent(in) :: dtime ! model time step (sec)
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'CalcIrrigNstepsPerDay'
    !-----------------------------------------------------------------------
    
    irrig_nsteps_per_day = ((this%params%irrig_length + (dtime - 1))/dtime)  ! round up

  end function CalcIrrigNstepsPerDay

  !-----------------------------------------------------------------------
  subroutine SetIrrigMethod(this, bounds)
    !
    ! !DESCRIPTION:
    ! Set this%irrig_method_patch based on surface dataset
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    !
    ! !LOCAL VARIABLES:
    integer :: p ! patch index
    integer :: g ! gridcell index
    integer :: m ! patch itype

    character(len=*), parameter :: subname = 'SetIrrigMethod'
    !-----------------------------------------------------------------------

    do p = bounds%begp,bounds%endp
       g = patch%gridcell(p)
       m = patch%itype(p)
       if (m >= lbound(irrig_method, 2) .and. m <= ubound(irrig_method, 2) &
            .and. pftcon%irrigated(m) == 1._r8) then
          this%irrig_method_patch(p) = irrig_method(g,m)
          ! ensure irrig_method is valid; if not set, use drip irrigation
          if(irrig_method(g,m) == irrig_method_unset) then
             this%irrig_method_patch(p) = this%params%irrig_method_default
          else if (irrig_method(g,m) /= irrig_method_drip .and. irrig_method(g,m) /= irrig_method_sprinkler) then
             write(iulog,*) subname //' invalid irrigation method specified'
             call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg='bad irrig_method '// &
                  errMsg(sourcefile, __LINE__))
          end if
       else
          this%irrig_method_patch(p) = this%params%irrig_method_default
       end if
    end do

  end subroutine SetIrrigMethod


  !-----------------------------------------------------------------------
  subroutine Restart(this, bounds, ncid, flag)
    !
    ! !DESCRIPTION:
    ! Handle restart of irrigation variables
    !
    ! !USES:
    use ncdio_pio        , only : file_desc_t, ncd_inqvdlen, ncd_double, ncd_int
    use restUtilMod
    !
    ! !ARGUMENTS:
    class(irrigation_type) :: this
    type(bounds_type), intent(in)    :: bounds 
    type(file_desc_t), intent(inout) :: ncid   ! netcdf id
    character(len=*) , intent(in)    :: flag   ! 'read', 'write' or 'define'
    !
    ! !LOCAL VARIABLES:
    logical :: do_io
    integer :: dimlen       ! dimension length
    integer :: nump_global  ! total number of patchs, globally
    integer :: err_code     ! error code
    logical :: readvar      ! determine if variable is on initial file

    character(len=*), parameter :: subname = 'Restart'
    !-----------------------------------------------------------------------
    
    ! Get expected total number of points, for later error checks
    call get_proc_global(np=nump_global)

    do_io = .true.
    readvar = .false.
    if (flag == 'read') then
       ! BACKWARDS_COMPATIBILITY
       ! On a read, confirm that this variable has the expected size; if not, don't read
       ! it (instead give it a default value). This is needed to support older initial
       ! conditions for which this variable had a different size.
       call ncd_inqvdlen(ncid, 'n_irrig_steps_left', 1, dimlen, err_code)
       if (dimlen /= nump_global) then
          do_io = .false.
       end if
    end if
    if (do_io) then
       call restartvar(ncid=ncid, flag=flag, varname='n_irrig_steps_left', xtype=ncd_int,  &
            dim1name='pft', &
            long_name='number of irrigation time steps left', units='#', &
            interpinic_flag='interp', readvar=readvar, data=this%n_irrig_steps_left_patch)
    end if
    if (flag=='read' .and. .not. readvar) then
       this%n_irrig_steps_left_patch = 0
    end if

    do_io = .true.
    readvar = .false.
    if (flag == 'read') then
       ! BACKWARDS_COMPATIBILITY
       ! On a read, confirm that this variable has the expected size; if not, don't read
       ! it (instead give it a default value). This is needed to support older initial
       ! conditions for which this variable had a different size.
       call ncd_inqvdlen(ncid, 'irrig_rate', 1, dimlen, err_code)
       if (dimlen /= nump_global) then
          do_io = .false.
       end if
    end if
    if (do_io) then
       call restartvar(ncid=ncid, flag=flag, varname='irrig_rate', xtype=ncd_double,  &
            dim1name='pft', &
            long_name='surface irrigation rate', units='mm/s', &
            interpinic_flag='interp', readvar=readvar, data=this%sfc_irrig_rate_patch)
    end if
    if (flag=='read' .and. .not. readvar) then
       this%sfc_irrig_rate_patch = 0.0_r8
    end if

    ! BACKWARDS_COMPATIBILITY(wjs, 2016-11-23) To support older restart files without an
    ! irrig_rate_demand field, get irrig_rate_demand from irrig_rate. I'm abusing the
    ! capability to specify multiple variable names here (even though irrig_rate isn't
    ! really an old version of irrig_rate_demand), rather than having code like 'if
    ! (.not. readvar) then ...', because the latter doesn't work when using init_interp.
    call restartvar(ncid=ncid, flag=flag, varname='irrig_rate_demand:irrig_rate', &
         xtype=ncd_double,  &
         dim1name='pft', &
         long_name='irrigation rate demand, neglecting surface water source limitation', &
         units='mm/s', &
         interpinic_flag='interp', readvar=readvar, data=this%irrig_rate_demand_patch)
  end subroutine Restart

  !-----------------------------------------------------------------------
  subroutine IrrigationClean(this)
    !
    ! !DESCRIPTION:
    ! Deallocate memory
    !
    ! !ARGUMENTS:
    class(irrigation_type), intent(inout) :: this
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'Clean'
    !-----------------------------------------------------------------------
    
    deallocate(this%qflx_irrig_demand_patch)
    deallocate(this%relsat_wilting_point_col)
    deallocate(this%relsat_target_col)
    deallocate(this%sfc_irrig_rate_patch)
    deallocate(this%irrig_rate_demand_patch)
    deallocate(this%irrig_method_patch)
    deallocate(this%n_irrig_steps_left_patch)

  end subroutine IrrigationClean


  ! ========================================================================
  ! Science routines
  ! ========================================================================
  
  !-----------------------------------------------------------------------
  subroutine CalcIrrigationFluxes(this, bounds, num_soilc, &
       filter_soilc, num_soilp, filter_soilp, &
       soilhydrology_inst, soilstate_inst, &
       water_inst)
    !
    ! !DESCRIPTION:
    ! Apply the irrigation computed by CalcIrrigationNeeded in order to set various fluxes
    !
    ! Sets irrigation withdrawal and application fluxes in the waterflux components of
    ! water_inst, for both bulk and tracers:
    ! - qflx_sfc_irrig_col
    ! - qflx_gw_uncon_irrig_lyr_col
    ! - qflx_gw_uncon_irrig_col
    ! - qflx_gw_con_irrig_col
    ! - qflx_irrig_drip_patch
    ! - qflx_irrig_sprinkler_patch
    !
    ! Should be called once, AND ONLY ONCE, per time step.
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds
    ! number of points in filter_soilc
    integer, intent(in) :: num_soilc
    ! column filter for soil
    integer, intent(in) :: filter_soilc(:)
    ! number of points in filter_soilp
    integer, intent(in) :: num_soilp
    ! patch filter for soil
    integer, intent(in) :: filter_soilp(:)
    type(soilhydrology_type) , intent(in)    :: soilhydrology_inst
    type(soilstate_type)     , intent(in)    :: soilstate_inst
    type(water_type)         , intent(inout) :: water_inst
    !
    ! !LOCAL VARIABLES:
    integer  :: p,fp  ! patch indices
    integer  :: c,fc  ! column indices
    integer  :: g     ! grid cell index
    integer  :: j     ! level index
    integer  :: i     ! tracer index
    real(r8) :: qflx_sfc_irrig_bulk_patch(bounds%begp:bounds%endp)      ! patch surface irrigation flux for bulk water (mm H2O/s)
    real(r8) :: qflx_gw_demand_bulk_patch(bounds%begp:bounds%endp) ! patch amount of irrigation groundwater demand for bulk water (mm H2O/s)
    real(r8) :: qflx_gw_demand_bulk_col(bounds%begc:bounds%endc) ! col amount of irrigation groundwater demand for bulk water (mm H2O/s)
    logical  :: is_bulk

    character(len=*), parameter :: subname = 'CalcIrrigationFluxes'

    !-----------------------------------------------------------------------

    ! This should be called exactly once per time step, so that the counter decrease
    ! works correctly.

    ! Note that these associated variables refer to the bulk instance
    associate( &
         begp => bounds%begp, &
         endp => bounds%endp, &
         begc => bounds%begc, &
         endc => bounds%endc  &
         )

    call this%CalcBulkWithdrawals(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
         soilhydrology_inst, soilstate_inst, water_inst%waterfluxbulk_inst, &
         qflx_sfc_irrig_bulk_patch = qflx_sfc_irrig_bulk_patch(begp:endp), &
         qflx_gw_demand_bulk_patch = qflx_gw_demand_bulk_patch(begp:endp), &
         qflx_gw_demand_bulk_col = qflx_gw_demand_bulk_col(begc:endc))

    do i = water_inst%tracers_beg, water_inst%tracers_end
       call this%CalcOneTracerWithdrawals(bounds, num_soilc, filter_soilc, &
            waterstatebulk_inst = water_inst%waterstatebulk_inst, &
            waterfluxbulk_inst = water_inst%waterfluxbulk_inst, &
            waterstate_tracer_inst = water_inst%bulk_and_tracers(i)%waterstate_inst, &
            waterflux_tracer_inst = water_inst%bulk_and_tracers(i)%waterflux_inst)
    end do

    if (this%params%use_groundwater_irrigation) then
       do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end
          call this%CalcTotalGWUnconIrrig(num_soilc, filter_soilc, &
               water_inst%bulk_and_tracers(i)%waterflux_inst)
       end do
    end if

    do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end
       is_bulk = (i == water_inst%i_bulk)
       call this%CalcApplicationFluxes(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
            waterflux_inst = water_inst%bulk_and_tracers(i)%waterflux_inst, &
            is_bulk = is_bulk, &
            qflx_sfc_irrig_bulk_patch = qflx_sfc_irrig_bulk_patch(begp:endp), &
            qflx_sfc_irrig_bulk_col = water_inst%waterfluxbulk_inst%qflx_sfc_irrig_col(begc:endc), &
            qflx_gw_demand_bulk_patch = qflx_gw_demand_bulk_patch(begp:endp), &
            qflx_gw_demand_bulk_col = qflx_gw_demand_bulk_col(begc:endc))
    end do

    end associate

  end subroutine CalcIrrigationFluxes

  !-----------------------------------------------------------------------
  subroutine CalcBulkWithdrawals(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
       soilhydrology_inst, soilstate_inst, waterfluxbulk_inst, &
       qflx_sfc_irrig_bulk_patch, qflx_gw_demand_bulk_patch, qflx_gw_demand_bulk_col)
    !
    ! !DESCRIPTION:
    ! Calculate irrigation withdrawals for bulk water
    !
    ! Sets / updates the following variables:
    ! - qflx_sfc_irrig_bulk_patch
    ! - qflx_gw_demand_bulk_patch
    ! - qflx_gw_demand_bulk_col
    ! - this%qflx_irrig_demand_patch
    ! - this%n_irrig_steps_left_patch
    ! - waterfluxbulk_inst%qflx_sfc_irrig_col
    ! - waterfluxbulk_inst%qflx_gw_uncon_irrig_lyr_col
    ! - waterfluxbulk_inst%qflx_gw_con_irrig_col
    !
    ! !ARGUMENTS:
    class(irrigation_type)   , intent(inout) :: this
    type(bounds_type)        , intent(in)    :: bounds
    integer                  , intent(in)    :: num_soilc       ! number of points in filter_soilc
    integer                  , intent(in)    :: filter_soilc(:) ! column filter for soil
    integer                  , intent(in)    :: num_soilp       ! number of points in filter_soilp
    integer                  , intent(in)    :: filter_soilp(:) ! patch filter for soil
    type(soilhydrology_type) , intent(in)    :: soilhydrology_inst
    type(soilstate_type)     , intent(in)    :: soilstate_inst
    type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst
    real(r8)                 , intent(inout) :: qflx_sfc_irrig_bulk_patch( bounds%begp: ) ! patch surface irrigation flux for bulk water (mm H2O/s)
    real(r8)                 , intent(inout) :: qflx_gw_demand_bulk_patch( bounds%begp: ) ! patch amount of irrigation groundwater demand for bulk water (mm H2O/s)
    real(r8)                 , intent(inout) :: qflx_gw_demand_bulk_col( bounds%begc: )   ! col amount of irrigation groundwater demand for bulk water (mm H2O/s)

    !
    ! !LOCAL VARIABLES:
    integer  :: p,fp  ! patch indices

    character(len=*), parameter :: subname = 'CalcBulkWithdrawals'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(qflx_sfc_irrig_bulk_patch) == [bounds%endp]), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(qflx_gw_demand_bulk_patch) == [bounds%endp]), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(qflx_gw_demand_bulk_col) == [bounds%endc]), sourcefile, __LINE__)

    associate( &
         begp => bounds%begp, &
         endp => bounds%endp, &
         begc => bounds%begc, &
         endc => bounds%endc, &

         qflx_sfc_irrig_col          => waterfluxbulk_inst%qflx_sfc_irrig_col          , & ! Output: [real(r8) (:)] col irrigation flux (mm H2O/s)
         qflx_gw_uncon_irrig_lyr_col => waterfluxbulk_inst%qflx_gw_uncon_irrig_lyr_col , & ! Output: [real(r8) (:,:) ] unconfined groundwater irrigation flux, separated by layer (mm H2O/s)
         qflx_gw_con_irrig_col       => waterfluxbulk_inst%qflx_gw_con_irrig_col         & ! Output: [real(r8) (:)] col confined groundwater irrigation flux (mm H2O/s)
         )

    do fp = 1, num_soilp
       p = filter_soilp(fp)
       
       if (this%n_irrig_steps_left_patch(p) > 0) then
          qflx_sfc_irrig_bulk_patch(p)     = this%sfc_irrig_rate_patch(p)
          this%qflx_irrig_demand_patch(p)  = this%irrig_rate_demand_patch(p)
          qflx_gw_demand_bulk_patch(p)     = &
               this%qflx_irrig_demand_patch(p) - qflx_sfc_irrig_bulk_patch(p)
          
          this%n_irrig_steps_left_patch(p) = this%n_irrig_steps_left_patch(p) - 1
       else
          qflx_sfc_irrig_bulk_patch(p)    = 0._r8
          this%qflx_irrig_demand_patch(p) = 0._r8
          qflx_gw_demand_bulk_patch(p)    = 0._r8
       end if
    end do

    call p2c (bounds, num_soilc, filter_soilc, &
         patcharr = qflx_sfc_irrig_bulk_patch(begp:endp), &
         colarr = qflx_sfc_irrig_col(begc:endc))

    call p2c (bounds, num_soilc, filter_soilc, &
         patcharr = qflx_gw_demand_bulk_patch(begp:endp), &
         colarr = qflx_gw_demand_bulk_col(begc:endc))

    if (this%params%use_groundwater_irrigation) then
       ! Supply as much of the remaining deficit as possible from groundwater irrigation
       call this%WrapCalcIrrigWithdrawals(bounds, num_soilc, filter_soilc, &
            soilhydrology_inst, soilstate_inst, &
            qflx_gw_demand = qflx_gw_demand_bulk_col(begc:endc), &
            qflx_gw_uncon_irrig_lyr = qflx_gw_uncon_irrig_lyr_col(begc:endc, :), &
            qflx_gw_con_irrig = qflx_gw_con_irrig_col(begc:endc))
    end if

    end associate

  end subroutine CalcBulkWithdrawals

  !-----------------------------------------------------------------------
  subroutine CalcOneTracerWithdrawals(this, bounds, num_soilc, filter_soilc, &
       waterstatebulk_inst, waterfluxbulk_inst, &
       waterstate_tracer_inst, waterflux_tracer_inst)
    !
    ! !DESCRIPTION:
    ! Calculate irrigation withdrawals for one water tracer
    !
    ! Sets the following variables:
    ! - waterflux_tracer_inst%qflx_sfc_irrig_col
    ! - waterflux_tracer_inst%qflx_gw_uncon_irrig_lyr_col
    ! - waterflux_tracer_inst%qflx_gw_con_irrig_col
    !
    ! !ARGUMENTS:
    class(irrigation_type)    , intent(in)    :: this
    type(bounds_type)         , intent(in)    :: bounds
    integer                   , intent(in)    :: num_soilc       ! number of points in filter_soilc
    integer                   , intent(in)    :: filter_soilc(:) ! column filter for soil
    type(waterstatebulk_type) , intent(in)    :: waterstatebulk_inst
    type(waterfluxbulk_type)  , intent(in)    :: waterfluxbulk_inst
    type(waterstate_type)     , intent(in)    :: waterstate_tracer_inst
    type(waterflux_type)      , intent(inout) :: waterflux_tracer_inst
    !
    ! !LOCAL VARIABLES:
    integer :: j  ! level index

    character(len=*), parameter :: subname = 'CalcOneTracerWithdrawals'
    !-----------------------------------------------------------------------

    associate( &
         begc => bounds%begc, &
         endc => bounds%endc &
         )

    ! BUG(wjs, 2018-12-17, ESCOMP/ctsm#512) Eventually, use actual source ratio for
    ! qflx_sfc_irrig_col rather than assuming a fixed ratio
    call CalcTracerFromBulkFixedRatio( &
         bulk   = waterfluxbulk_inst%qflx_sfc_irrig_col(begc:endc), &
         ratio  = waterflux_tracer_inst%info%get_ratio(), &
         tracer = waterflux_tracer_inst%qflx_sfc_irrig_col(begc:endc))

    if (this%params%use_groundwater_irrigation) then

       call CalcTracerFromBulk( &
            subgrid_level = subgrid_level_column, &
            lb            = begc, &
            num_pts       = num_soilc, &
            filter_pts    = filter_soilc, &
            bulk_source   = waterstatebulk_inst%wa_col(begc:endc), &
            bulk_val      = waterfluxbulk_inst%qflx_gw_con_irrig_col(begc:endc), &
            tracer_source = waterstate_tracer_inst%wa_col(begc:endc), &
            tracer_val    = waterflux_tracer_inst%qflx_gw_con_irrig_col(begc:endc))
       do j = 1, nlevsoi
          call CalcTracerFromBulk( &
               subgrid_level = subgrid_level_column, &
               lb            = begc, &
               num_pts       = num_soilc, &
               filter_pts    = filter_soilc, &
               bulk_source   = waterstatebulk_inst%h2osoi_liq_col(begc:endc, j), &
               bulk_val      = waterfluxbulk_inst%qflx_gw_uncon_irrig_lyr_col(begc:endc, j), &
               tracer_source = waterstate_tracer_inst%h2osoi_liq_col(begc:endc, j), &
               tracer_val    = waterflux_tracer_inst%qflx_gw_uncon_irrig_lyr_col(begc:endc, j))
       end do

    end if

    end associate

  end subroutine CalcOneTracerWithdrawals

  !-----------------------------------------------------------------------
  subroutine CalcTotalGWUnconIrrig(this, num_soilc, filter_soilc, &
       waterflux_inst)
    !
    ! !DESCRIPTION:
    ! Calculate total irrigation withdrawal flux from the unconfined aquifer, for either bulk
    ! or one water tracer
    !
    ! Sets the following variables:
    ! - waterflux_inst%qflx_gw_uncon_irrig_col
    !
    ! !ARGUMENTS:
    class(irrigation_type)    , intent(in)    :: this
    integer                   , intent(in)    :: num_soilc       ! number of points in filter_soilc
    integer                   , intent(in)    :: filter_soilc(:) ! column filter for soil
    class(waterflux_type)     , intent(inout) :: waterflux_inst
    !
    ! !LOCAL VARIABLES:
    integer :: c, fc  ! column indices
    integer :: j      ! level index

    character(len=*), parameter :: subname = 'CalcTotalGWUnconIrrig'
    !-----------------------------------------------------------------------

    do fc = 1, num_soilc
       c = filter_soilc(fc)
       waterflux_inst%qflx_gw_uncon_irrig_col(c) = 0._r8
    end do
    do j = 1, nlevsoi
       do fc = 1, num_soilc
          c = filter_soilc(fc)
          waterflux_inst%qflx_gw_uncon_irrig_col(c) = &
               waterflux_inst%qflx_gw_uncon_irrig_col(c) + &
               waterflux_inst%qflx_gw_uncon_irrig_lyr_col(c,j)
       end do
    end do

  end subroutine CalcTotalGWUnconIrrig

  !-----------------------------------------------------------------------
  subroutine CalcApplicationFluxes(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, &
       waterflux_inst, is_bulk, &
       qflx_sfc_irrig_bulk_patch, qflx_sfc_irrig_bulk_col, &
       qflx_gw_demand_bulk_patch, qflx_gw_demand_bulk_col)
    !
    ! !DESCRIPTION:
    ! Calculate irrigation application fluxes for either bulk water or a single tracer
    !
    ! Sets:
    ! - waterflux_inst%qflx_irrig_drip_patch
    ! - waterflux_inst%qflx_irrig_sprinkler_patch
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(in)    :: this
    type(bounds_type)      , intent(in)    :: bounds
    integer                , intent(in)    :: num_soilc       ! number of points in filter_soilc
    integer                , intent(in)    :: filter_soilc(:) ! column filter for soil
    integer                , intent(in)    :: num_soilp       ! number of points in filter_soilp
    integer                , intent(in)    :: filter_soilp(:) ! patch filter for soil
    class(waterflux_type)  , intent(inout) :: waterflux_inst
    logical                , intent(in)    :: is_bulk         ! whether this call is being made for bulk water (is_bulk true) or a tracer (is_bulk false)
    real(r8)               , intent(in)    :: qflx_sfc_irrig_bulk_patch( bounds%begp: ) ! patch surface irrigation flux for bulk water (mm H2O/s)
    real(r8)               , intent(in)    :: qflx_sfc_irrig_bulk_col( bounds%begc: )   ! col surface irrigation flux for bulk water (mm H2O/s)
    real(r8)               , intent(in)    :: qflx_gw_demand_bulk_patch( bounds%begp: ) ! patch amount of irrigation groundwater demand for bulk water (mm H2O/s)
    real(r8)               , intent(in)    :: qflx_gw_demand_bulk_col( bounds%begc: )   ! col amount of irrigation groundwater demand for bulk water (mm H2O/s)
    !
    ! !LOCAL VARIABLES:
    integer  :: fc, c  ! column indices
    integer  :: fp, p  ! patch indices
    real(r8) :: qflx_gw_irrig_withdrawn_col(bounds%begc:bounds%endc) ! col total amount of irrigation withdrawn from groundwater (mm H2O/s)
    real(r8) :: qflx_sfc_irrig        ! amount of surface irrigation for a single patch (mm H2O/s)
    real(r8) :: qflx_gw_irrig_applied ! amount of groundwater irrigation for a single patch (mm H2O/s)
    real(r8) :: qflx_irrig_tot        ! total amount of irrigation for a single patch (mm H2O/s)

    character(len=*), parameter :: subname = 'CalcApplicationFluxes'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(qflx_sfc_irrig_bulk_patch) == [bounds%endp]), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(qflx_sfc_irrig_bulk_col) == [bounds%endc]), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(qflx_gw_demand_bulk_patch) == [bounds%endp]), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(qflx_gw_demand_bulk_col) == [bounds%endc]), sourcefile, __LINE__)

    do fc = 1, num_soilc
       c = filter_soilc(fc)
       ! Note that qflx_gw_uncon_irrig_col and qflx_gw_con_irrig_col will remain at
       ! their cold start values of 0 if use_groundwater_irrigation is false, so this
       ! sum is safe to do in all cases.
       qflx_gw_irrig_withdrawn_col(c) = &
            waterflux_inst%qflx_gw_uncon_irrig_col(c) + &
            waterflux_inst%qflx_gw_con_irrig_col(c)
    end do

    do fp = 1, num_soilp
       p = filter_soilp(fp)
       c = patch%column(p)

       if (is_bulk) then
          qflx_sfc_irrig = qflx_sfc_irrig_bulk_patch(p)
       else
          if (qflx_sfc_irrig_bulk_col(c) > 0._r8) then
             ! We have col-level sfc irrig for each water tracer, but only have
             ! patch-level for the bulk. To get the patch-level tracer flux, we need
             ! to scale the column-level flux by the ratio of patch-to-column-level
             ! flux for the bulk.
             !
             ! Note that this is similar to the scaling that we do for groundwater
             ! irrigation below.
             qflx_sfc_irrig = waterflux_inst%qflx_sfc_irrig_col(c) * &
                  (qflx_sfc_irrig_bulk_patch(p) / qflx_sfc_irrig_bulk_col(c))
          else
             if (qflx_sfc_irrig_bulk_patch(p) > 0._r8 .or. &
                  waterflux_inst%qflx_sfc_irrig_col(c) > 0._r8) then
                write(iulog,*) 'If qflx_sfc_irrig_bulk_col <= 0, ' // &
                     'expect qflx_sfc_irrig_bulk_patch = waterflux_inst%qflx_sfc_irrig_col = 0'
                write(iulog,*) 'qflx_sfc_irrig_bulk_col = ', qflx_sfc_irrig_bulk_col(c)
                write(iulog,*) 'qflx_sfc_irrig_bulk_patch = ', qflx_sfc_irrig_bulk_patch(p)
                write(iulog,*) 'waterflux_inst%qflx_sfc_irrig_col = ', &
                     waterflux_inst%qflx_sfc_irrig_col(c)
                call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, &
                     msg = 'If qflx_sfc_irrig_bulk_col <= 0, ' // &
                     'expect qflx_sfc_irrig_bulk_patch = waterflux_inst%qflx_sfc_irrig_col = 0', &
                     additional_msg = errMsg(sourcefile, __LINE__))
             end if

             qflx_sfc_irrig = 0._r8
          end if
       end if

       if (qflx_gw_demand_bulk_col(c) > 0._r8) then
          ! Distribute the withdrawn groundwater irrigation to each patch according to
          ! its demand relative to the total column demand.
          !
          ! Note that this expression could be reformulated to:
          !   qflx_gw_irrig_applied = qflx_gw_demand_bulk_patch(p) * &
          !        (qflx_gw_irrig_withdrawn_col(c) / qflx_gw_demand_bulk_col(c))
          !
          ! It may be more intuitive that the above is correct, at least for bulk water:
          ! this shows that we're down-weighting each patch's demand evenly, based on the
          ! ration of withdrawn water to demand in the given column.
          !
          ! Note that this is similar to the scaling that we do for surface irrigation
          ! above.
          qflx_gw_irrig_applied = qflx_gw_irrig_withdrawn_col(c) * &
               (qflx_gw_demand_bulk_patch(p) / qflx_gw_demand_bulk_col(c))
       else
          if (qflx_gw_demand_bulk_patch(p) > 0._r8 .or. &
               qflx_gw_irrig_withdrawn_col(c) > 0._r8) then
             write(iulog,*) 'If qflx_gw_demand_bulk_col <= 0, expect qflx_gw_demand_bulk_patch = qflx_gw_irrig_withdrawn_col = 0'
             write(iulog,*) 'qflx_gw_demand_bulk_col = ', qflx_gw_demand_bulk_col(c)
             write(iulog,*) 'qflx_gw_demand_bulk_patch = ', qflx_gw_demand_bulk_patch(p)
             write(iulog,*) 'qflx_gw_irrig_withdrawn_col = ', qflx_gw_irrig_withdrawn_col(c)
             call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, &
                  msg = 'If qflx_gw_demand_bulk_col <= 0, expect qflx_gw_demand_bulk_patch = qflx_gw_irrig_withdrawn_col = 0', &
                  additional_msg = errMsg(sourcefile, __LINE__))
          end if

          qflx_gw_irrig_applied = 0._r8
       end if

       qflx_irrig_tot = qflx_sfc_irrig + qflx_gw_irrig_applied

       ! Set drip/sprinkler irrigation based on irrigation method from input data
       waterflux_inst%qflx_irrig_drip_patch(p)      = 0._r8
       waterflux_inst%qflx_irrig_sprinkler_patch(p) = 0._r8

       if(this%irrig_method_patch(p) == irrig_method_drip) then
          waterflux_inst%qflx_irrig_drip_patch(p)      = qflx_irrig_tot
       else if(this%irrig_method_patch(p) == irrig_method_sprinkler) then
          waterflux_inst%qflx_irrig_sprinkler_patch(p) = qflx_irrig_tot
       else
          call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, &
               msg=' ERROR: irrig_method_patch set to invalid value ' // &
               errMsg(sourcefile, __LINE__))
       endif

    end do

  end subroutine CalcApplicationFluxes

  !-----------------------------------------------------------------------
  subroutine WrapCalcIrrigWithdrawals(this, bounds, num_soilc, filter_soilc, &
       soilhydrology_inst, soilstate_inst, &
       qflx_gw_demand, &
       qflx_gw_uncon_irrig_lyr, &
       qflx_gw_con_irrig)
    !
    ! !DESCRIPTION:
    ! Wrap the call to the external subroutine, CalcIrrigWithdrawals
    !
    ! The purpose of this wrapping is to allow unit tests to override the behavior
    !
    ! !ARGUMENTS:
    class(irrigation_type)   , intent(in)    :: this
    type(bounds_type)        , intent(in)    :: bounds
    integer                  , intent(in)    :: num_soilc       ! number of points in filter_soilc
    integer                  , intent(in)    :: filter_soilc(:) ! column filter for soil
    type(soilhydrology_type) , intent(in)    :: soilhydrology_inst
    type(soilstate_type)     , intent(in)    :: soilstate_inst

    real(r8) , intent(in)    :: qflx_gw_demand( bounds%begc: )              ! groundwater irrigation demand (mm H2O/s)
    real(r8) , intent(inout) :: qflx_gw_uncon_irrig_lyr( bounds%begc:, 1: ) ! unconfined aquifer groundwater irrigation withdrawal flux, separated by layer (mm H2O/s)
    real(r8) , intent(inout) :: qflx_gw_con_irrig( bounds%begc: )           ! total confined aquifer groundwater irrigation withdrawal flux (mm H2O/s)
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'WrapCalcIrrigWithdrawals'
    !-----------------------------------------------------------------------

    call CalcIrrigWithdrawals( &
         bounds = bounds, &
         num_soilc = num_soilc, &
         filter_soilc = filter_soilc, &
         soilhydrology_inst = soilhydrology_inst, &
         soilstate_inst = soilstate_inst, &
         qflx_gw_demand = qflx_gw_demand(bounds%begc:bounds%endc), &
         qflx_gw_uncon_irrig_lyr = qflx_gw_uncon_irrig_lyr(bounds%begc:bounds%endc, :), &
         qflx_gw_con_irrig = qflx_gw_con_irrig(bounds%begc:bounds%endc))

  end subroutine WrapCalcIrrigWithdrawals

  !-----------------------------------------------------------------------
  subroutine CalcIrrigationNeeded(this, bounds, num_exposedvegp, filter_exposedvegp, &
       elai, t_soisno, eff_porosity, h2osoi_liq, volr, rof_prognostic)
    !
    ! !DESCRIPTION:
    ! Calculate whether and how much irrigation is needed for each column. However, this
    ! does NOT actually set the irrigation flux.
    !
    ! !USES:
    use shr_const_mod      , only : SHR_CONST_TKFRZ
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(inout) :: this
    type(bounds_type)      , intent(in)    :: bounds

    ! number of points in filter_exposedvegp
    integer, intent(in) :: num_exposedvegp

    ! patch filter for non-snow-covered veg
    integer, intent(in) :: filter_exposedvegp(:)

    ! one-sided leaf area index with burying by snow [patch]
    real(r8), intent(in) :: elai( bounds%begp: )

    ! col soil temperature (K) [col, nlevgrnd] (note that this does NOT contain the snow levels)
    real(r8), intent(in) :: t_soisno( bounds%begc: , 1: )

    ! effective porosity (0 to 1) [col, nlevgrnd]
    real(r8), intent(in) :: eff_porosity( bounds%begc: , 1: )

    ! column liquid water (kg/m2) [col, nlevgrnd] (note that this does NOT contain the snow levels)
    real(r8), intent(in) :: h2osoi_liq( bounds%begc: , 1: )

    ! river water volume (m3) (ignored if rof_prognostic is .false.)
    real(r8), intent(in) :: volr( bounds%begg: )

    ! whether we're running with a prognostic ROF component; this is needed to determine
    ! whether we can limit irrigation based on river volume.
    logical, intent(in) :: rof_prognostic

    !
    ! !LOCAL VARIABLES:
    integer :: fp   ! patch filter index
    integer :: fc   ! column filter index
    integer :: p    ! patch index
    integer :: c    ! column index
    integer :: g    ! gridcell index
    integer :: j    ! level

    ! Filter for columns where we need to check for irrigation
    type(filter_col_type) :: check_for_irrig_col_filter

    ! target soil moisture for this layer [kg/m2]
    real(r8) :: h2osoi_liq_target

    ! soil moisture at wilting point for this layer [kg/m2]
    real(r8) :: h2osoi_liq_wilting_point

    ! Total of h2osoi down to the depth of irrigation in each column [kg/m2]
    real(r8) :: h2osoi_liq_tot(bounds%begc:bounds%endc)

    ! Total of h2osoi_liq_target down to the depth of irrigation in each column [kg/m2]
    real(r8) :: h2osoi_liq_target_tot(bounds%begc:bounds%endc)

    ! Total of h2osoi_liq at wilting point down to the depth of irrigation in each column
    ! [kg/m2]
    real(r8) :: h2osoi_liq_wilting_point_tot(bounds%begc:bounds%endc)

    ! h2osoi_liq at the threshold for irrigation in this column [kg/m2]
    real(r8) :: h2osoi_liq_at_threshold

    ! difference between desired soil moisture level for each column and current soil
    ! moisture level [kg/m2] [i.e., mm]
    real(r8) :: deficit(bounds%begc:bounds%endc)

    ! deficit limited by river volume [kg/m2] [i.e., mm]
    real(r8) :: deficit_volr_limited(bounds%begc:bounds%endc)

    ! where do we need to check soil moisture to see if we need to irrigate?
    logical  :: check_for_irrig_patch(bounds%begp:bounds%endp)
    logical  :: check_for_irrig_col(bounds%begc:bounds%endc)

    ! set to true once we have reached the max allowable depth for irrigation in a given
    ! column
    logical  :: reached_max_depth(bounds%begc:bounds%endc)

    ! Whether we should limit deficits by available volr
    logical :: limit_irrigation

    character(len=*), parameter :: subname = 'CalcIrrigationNeeded'
    !-----------------------------------------------------------------------
    
    ! Enforce expected array sizes
    SHR_ASSERT_ALL_FL((ubound(elai) == (/bounds%endp/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(t_soisno) == (/bounds%endc, nlevgrnd/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(eff_porosity) == (/bounds%endc, nlevgrnd/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(h2osoi_liq) == (/bounds%endc, nlevgrnd/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(volr) == (/bounds%endg/)), sourcefile, __LINE__)


    ! Determine if irrigation is needed (over irrigated soil columns)
    
    ! First, determine in what grid cells we need to bother 'measuring' soil water, to see
    ! if we need irrigation
    check_for_irrig_col(bounds%begc:bounds%endc) = .false.
    do fp = 1, num_exposedvegp
       p = filter_exposedvegp(fp)
       g = patch%gridcell(p)

       check_for_irrig_patch(p) = this%PointNeedsCheckForIrrig( &
            pft_type=patch%itype(p), elai=elai(p), &
            londeg=grc%londeg(g))
       if (check_for_irrig_patch(p)) then
          c = patch%column(p)
          check_for_irrig_col(c) = .true.
       end if
    end do

    check_for_irrig_col_filter = col_filter_from_logical_array(bounds, &
         check_for_irrig_col(bounds%begc:bounds%endc))

    ! Initialize some variables
    do fc = 1, check_for_irrig_col_filter%num
       c = check_for_irrig_col_filter%indices(fc)

       reached_max_depth(c) = .false.
       h2osoi_liq_tot(c) = 0._r8
       h2osoi_liq_target_tot(c) = 0._r8
       h2osoi_liq_wilting_point_tot(c) = 0._r8
    end do

    ! Now 'measure' soil water for the grid cells identified above and see if the soil is
    ! dry enough to warrant irrigation
    do j = 1,nlevsoi
       do fc = 1, check_for_irrig_col_filter%num
          c = check_for_irrig_col_filter%indices(fc)

          if (.not. reached_max_depth(c)) then
             if (col%z(c,j) > this%params%irrig_depth) then
                reached_max_depth(c) = .true.
             else if (j > col%nbedrock(c)) then
                reached_max_depth(c) = .true.
             else if (t_soisno(c,j) <= SHR_CONST_TKFRZ) then
                ! if level L was frozen, then we don't look at any levels below L
                reached_max_depth(c) = .true.
             else
                h2osoi_liq_tot(c) = h2osoi_liq_tot(c) + h2osoi_liq(c,j)

                h2osoi_liq_target = this%RelsatToH2osoi( &
                     relsat = this%relsat_target_col(c,j), &
                     eff_porosity = eff_porosity(c,j), &
                     dz = col%dz(c,j))
                h2osoi_liq_target_tot(c) = h2osoi_liq_target_tot(c) + &
                     h2osoi_liq_target

                h2osoi_liq_wilting_point = this%RelsatToH2osoi( &
                     relsat = this%relsat_wilting_point_col(c,j), &
                     eff_porosity = eff_porosity(c,j), &
                     dz = col%dz(c,j))
                h2osoi_liq_wilting_point_tot(c) = h2osoi_liq_wilting_point_tot(c) + &
                     h2osoi_liq_wilting_point
             end if
          end if     ! if (.not. reached_max_depth(c))
       end do        ! do fc
    end do           ! do j

    ! Compute deficits
    ! First initialize deficits to 0 everywhere; this is needed for later averaging up to gridcell
    deficit(bounds%begc:bounds%endc) = 0._r8
    do fc = 1, check_for_irrig_col_filter%num
       c = check_for_irrig_col_filter%indices(fc)

       h2osoi_liq_at_threshold = h2osoi_liq_wilting_point_tot(c) + &
            this%params%irrig_threshold_fraction * &
            (h2osoi_liq_target_tot(c) - h2osoi_liq_wilting_point_tot(c))
       if (h2osoi_liq_tot(c) < h2osoi_liq_at_threshold) then
          deficit(c) = h2osoi_liq_target_tot(c) - h2osoi_liq_tot(c)
          ! deficit shouldn't be less than 0: if it is, that implies that the
          ! irrigation target is less than the irrigation threshold, which is not
          ! supposed to happen
          if (deficit(c) < 0._r8) then
             write(iulog,*) subname//' ERROR: deficit < 0'
             write(iulog,*) 'This implies that irrigation target is less than irrigatio&
                  &n threshold, which should never happen'
             call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg='deficit < 0 '// &
                  errMsg(sourcefile, __LINE__))
          end if
       else
          ! We're above the threshold - so don't irrigate
          deficit(c) = 0._r8
       end if
    end do

    ! Limit deficits by available volr, if desired. Note that we cannot do this limiting
    ! if running without a prognostic river model, since we need river volume for this
    ! limiting.
    !
    ! NOTE(wjs, 2016-11-22) In principle we could base this on rof_present rather than
    ! rof_prognostic, but that would depend on the data runoff (drof) model sending river
    ! volume, which it currently does not.
    limit_irrigation = (this%params%limit_irrigation_if_rof_enabled .and. rof_prognostic)
    if (limit_irrigation) then
       call this%CalcDeficitVolrLimited( &
            bounds = bounds, &
            check_for_irrig_col_filter = check_for_irrig_col_filter, &
            deficit = deficit(bounds%begc:bounds%endc), &
            volr = volr(bounds%begg:bounds%endg), &
            deficit_volr_limited = deficit_volr_limited(bounds%begc:bounds%endc))
    else
       deficit_volr_limited(bounds%begc:bounds%endc) = deficit(bounds%begc:bounds%endc)
    end if

    ! Convert deficits to irrigation rate
    do fp = 1, num_exposedvegp
       p = filter_exposedvegp(fp)
       c = patch%column(p)

       if (check_for_irrig_patch(p)) then

          ! Convert units from mm to mm/sec
          this%sfc_irrig_rate_patch(p) = deficit_volr_limited(c) / &
               (this%dtime*this%irrig_nsteps_per_day)
          this%irrig_rate_demand_patch(p) = deficit(c) / &
               (this%dtime*this%irrig_nsteps_per_day)

          ! n_irrig_steps_left(p) > 0 is ok even if irrig_rate(p) ends up = 0
          ! in this case, we'll irrigate by 0 for the given number of time steps
          this%n_irrig_steps_left_patch(p) = this%irrig_nsteps_per_day
       end if
    end do

  end subroutine CalcIrrigationNeeded

  !-----------------------------------------------------------------------
  function PointNeedsCheckForIrrig(this, pft_type, elai, londeg) &
       result(check_for_irrig)
    !
    ! !DESCRIPTION:
    ! Determine whether a given patch needs to be checked for irrigation now.
    !
    ! !USES:
    use clm_time_manager, only : get_local_time
    use pftconMod       , only : pftcon
    !
    ! !ARGUMENTS:
    logical :: check_for_irrig  ! function result
    class(irrigation_type), intent(in) :: this
    integer , intent(in) :: pft_type  ! type of pft in this patch
    real(r8), intent(in) :: elai      ! one-sided leaf area index with burying by snow
    real(r8), intent(in) :: londeg    ! longitude (degrees)
    !
    ! !LOCAL VARIABLES:
    ! number of seconds since the prescribed irrigation start time
    integer  :: seconds_since_irrig_start_time

    character(len=*), parameter :: subname = 'PointNeedsCheckForIrrig'
    !-----------------------------------------------------------------------
    
    if (pftcon%irrigated(pft_type) == 1._r8 .and. &
         elai > this%params%irrig_min_lai) then
       ! see if it's the right time of day to start irrigating:
       seconds_since_irrig_start_time = get_local_time( londeg, starttime=this%params%irrig_start_time, offset=-this%dtime )
       if (seconds_since_irrig_start_time < this%dtime) then
          check_for_irrig         = .true.
       else
          check_for_irrig    = .false.
       end if
    else
       check_for_irrig       = .false.
    end if

  end function PointNeedsCheckForIrrig

  !-----------------------------------------------------------------------
  subroutine CalcDeficitVolrLimited(this, bounds, check_for_irrig_col_filter, &
       deficit, volr, deficit_volr_limited)
    !
    ! !DESCRIPTION:
    ! Calculates deficit limited by river volume for each column.
    !
    ! The output array (deficit_volr_limited) is set to 0 outside the given filter.
    !
    ! If the deficit is lower than the volr-based threshold for a given gridcell, then the
    ! volr-limited deficit is simply the original deficit; if the deficit is higher than
    ! the volr-based threshold for a gridcell, then the volr-limited deficit for each
    ! column in the gridcell is equal to:
    !    (original deficit in the column) * (volr-based threshold)/(gridcell-level deficit)
    !
    ! The logic here relies on the fact that all irrigated columns in a given grid cell
    ! will be irrigated at the same time of day. (As opposed to, say, some being irrigated
    ! at 4 AM and others being irrigated at 6 AM: in that case we couldn't just average
    ! this time step's column-level deficits to the gridcell to get the total daily
    ! irrigation deficit: we'd need to account for the fact that we have some irrigation
    ! commitment from earlier deficits, etc.)
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    class(irrigation_type) , intent(in) :: this
    type(bounds_type)      , intent(in) :: bounds

    ! filter of columns where we need to check for irrigation
    type(filter_col_type)  , intent(in) :: check_for_irrig_col_filter

    ! original deficit: the amount by which we would irrigate if there were no volr
    ! limitation [kg/m2] [i.e., mm]
    !
    ! Should be set for all columns, even those outside the filter, so that averaging to
    ! grid cell happens properly
    real(r8), intent(in) :: deficit( bounds%begc: )

    ! river water volume [m3]
    real(r8), intent(in) :: volr( bounds%begg: )

    ! deficit limited by river volume [kg/m2] [i.e., mm]
    real(r8), intent(out) :: deficit_volr_limited( bounds%begc: )
    !
    ! !LOCAL VARIABLES:
    integer :: fc ! column filter index
    integer :: c  ! column index
    integer :: g  ! gridcell index

    ! Deficit averaged up to gridcell level [kg/m2] [i.e., mm]
    real(r8) :: deficit_grc(bounds%begg:bounds%endg)

    real(r8) :: available_volr ! volr available for withdrawal [m3]

    real(r8) :: max_deficit_supported_by_volr ! [kg/m2] [i.e., mm]

    ! ratio of deficit_volr_limited to deficit for each grid cell
    real(r8) :: deficit_limited_ratio_grc(bounds%begg:bounds%endg)

    character(len=*), parameter :: subname = 'CalcDeficitVolrLimited'
    !-----------------------------------------------------------------------

    SHR_ASSERT_ALL_FL((ubound(deficit) == (/bounds%endc/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(volr) == (/bounds%endg/)), sourcefile, __LINE__)
    SHR_ASSERT_ALL_FL((ubound(deficit_volr_limited) == (/bounds%endc/)), sourcefile, __LINE__)

    call c2g(bounds, &
         carr = deficit(bounds%begc:bounds%endc), &
         garr = deficit_grc(bounds%begg:bounds%endg), &
         c2l_scale_type = 'unity', &
         l2g_scale_type = 'unity')

    do g = bounds%begg, bounds%endg
       if (volr(g) > 0._r8) then
          available_volr = volr(g) * (1._r8 - this%params%irrig_river_volume_threshold)
          max_deficit_supported_by_volr = available_volr / grc%area(g) * m3_over_km2_to_mm
       else
          ! Ensure that negative volr is treated the same as 0 volr
          max_deficit_supported_by_volr = 0._r8
       end if

       if (deficit_grc(g) > max_deficit_supported_by_volr) then
          ! inadequate river storage, adjust irrigation demand
          deficit_limited_ratio_grc(g) = max_deficit_supported_by_volr / deficit_grc(g)
       else
          ! adequate river storage, no adjustment to irrigation demand
          deficit_limited_ratio_grc(g) = 1._r8
       end if
    end do

    deficit_volr_limited(bounds%begc:bounds%endc) = 0._r8
    do fc = 1, check_for_irrig_col_filter%num
       c = check_for_irrig_col_filter%indices(fc)

       g = col%gridcell(c)
       deficit_volr_limited(c) = deficit(c) * deficit_limited_ratio_grc(g)
    end do

  end subroutine CalcDeficitVolrLimited


  !-----------------------------------------------------------------------
  pure function RelsatToH2osoi(this, relsat, eff_porosity, dz) result(h2osoi_liq)
    !
    ! !DESCRIPTION:
    ! For a given column and layer, convert relative saturation to kg/m2 water
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8) :: h2osoi_liq   ! function result (kg/m2)
    class(irrigation_type), intent(in) :: this
    real(r8), intent(in) :: relsat       ! relative saturation (0 to 1)
    real(r8), intent(in) :: eff_porosity ! effective porosity (0 to 1)
    real(r8), intent(in) :: dz           ! level thickness (m)
    !
    ! !LOCAL VARIABLES:
    real(r8) :: vol_liq  ! partial volume of liquid water in layer [0 to 1]

    character(len=*), parameter :: subname = 'RelsatToH2osoi'
    !-----------------------------------------------------------------------

    vol_liq    = eff_porosity * relsat
    h2osoi_liq = vol_liq * denh2o * dz

  end function RelsatToH2osoi

  !-----------------------------------------------------------------------
  function UseGroundwaterIrrigation(this)
    !
    ! !DESCRIPTION:
    ! Returns true if we're using groundwater irrigation in this run
    !
    ! !ARGUMENTS:
    implicit none
    class(irrigation_type), intent(in) :: this

    logical :: UseGroundwaterIrrigation  ! function result
    !-----------------------------------------------------------------------

    UseGroundwaterIrrigation = this%params%use_groundwater_irrigation
    
  end function UseGroundwaterIrrigation
  
end module IrrigationMod
